#include "Random.hpp"

unsigned int Random::R250Index;
unsigned int Random::R250Buffer[250];
unsigned int Random::InitialSeed = Random::MAX_RANDOM;
unsigned int Random::CurrentSeed = Random::MAX_RANDOM;
bool Random::GaussianComputed = false;
double Random::LastGaussian = 0.0;

Random::Random()
{
    initialize( time( NULL ) );
}

Random::Random( unsigned int newSeed )
{
    initialize( newSeed );
}

Random::Random( const Random &other )
{
    initialize( other.initialSeed() );
}

Random::~Random()
{

}

Random &Random::operator=( const Random &other )
                         {
    initialize( other.initialSeed() );
    return *this;
}

void Random::initialize( unsigned int newSeed )
{
    InitialSeed = newSeed;
    CurrentSeed = newSeed;
    
    R250Index = 0;

    for( int j = 0; j < 250; j++ ) {     /* Fill the R250 buffer with 15-bit values */
        R250Buffer[j] = nextLCPNG();
    }

    for( unsigned int j = 0; j < 250; j++ ) {    /* Set some of the MS bits to 1 */
        if( nextLCPNG() > 16384) {
            R250Buffer[j] |= 0x8000;
        }
    }

    unsigned int msb = 0x8000;                   /* To turn on the diagonal bit   */
    unsigned int mask = 0xFFFF;                  /* To turn off the leftmost bits */
    for( unsigned int j = 0; j < 16; j++ ) {
        unsigned int k = 11 * j + 3;             /* Select a word to operate on        */
        R250Buffer[k] &= mask;             /* Turn off bits left of the diagonal */
        R250Buffer[k] |= msb;              /* Turn on the diagonal bit           */
        mask >>= 1;
        msb >>= 1;
    }
}

/* return uniformly distributed random on [0, MAX_RANDOM] */
unsigned int Random::next()
{
    int j;
    int new_rand;

    if( R250Index >= 147 ) {
        j = R250Index - 147;     /* Wrap pointer around */
    }
    else {
        j = R250Index + 103;
    }

    new_rand = R250Buffer[R250Index] ^= R250Buffer[j];

    if( R250Index >= 249 ) {      /* Increment pointer for next time */
        R250Index = 0;
    }
    else {
        R250Index++;
    }

    return new_rand;   /* Return a number */
}

/* return uniformly distributed random on [0, max] */
unsigned int Random::next( unsigned int max )
{
    return (int)( doubleNext() * ( max + 1 ) );
}

/* return uniformly distributed random on [min, max] */
int Random::next( int min, int max )
{
    return next( max - min ) + min;
}

/* return uniformly distributed random on [0.0, 1.0] */
double Random::doubleNext()
{
    return (double)next() / (double)MAX_RANDOM;
}

/* return uniformly distributed random on [0.0, max] */
double Random::doubleNext( double max )
{
    return doubleNext() * max;
}

/* return uniformly distributed random on [min, max] */
double Random::doubleNext( double min, double max )
{
    return doubleNext( max - min ) + min;
}

/* return gaussian distributed random for mean:0.0 and standard deviation: 1.0 */
double Random::doubleNextGaussian()
{
    // Polar form of Box-Muller transformation
    // http://www.taygeta.com/random/gaussian.html
    if( GaussianComputed ) {
        GaussianComputed = false;
        return LastGaussian;
    }
    else {
        double x1, x2, w;

        do {
            x1 = ( 2.0 * doubleNext() ) - 1.0;
            x2 = ( 2.0 * doubleNext() ) - 1.0;
            w = ( x1 * x1 ) + ( x2 * x2 );
        } while ( w >= 1.0 );

        w = sqrt( (-2.0 * log( w ) ) / w );

        GaussianComputed = true;
        LastGaussian = x1 * w;

        return x2 * w;
    }
}

/* return gaussian distributed random for mean:0.0 and standard deviation:standardDeviation  */
double Random::doubleNextGaussian( double standardDeviation )
{
    return doubleNextGaussian() * standardDeviation;
}

/* return gaussian distributed random for mean:mean and standard deviation:standardDeviation */
double Random::doubleNextGaussian( double mean, double standardDeviation )
{
    return doubleNextGaussian( standardDeviation ) + mean;
}

/* return gaussian distributed random for mean:0 and standard deviation:1 */
int Random::nextGaussian()
{
    return (unsigned int)( doubleNextGaussian() );
}

/* return gaussian distributed random for mean:0 and standard deviation:standardDeviation */
int Random::nextGaussian( int standardDeviation )
{
    return (unsigned int)(doubleNextGaussian() * standardDeviation);
}

/* return gaussian distributed random for mean:mean and standard deviation:standardDeviation */
int Random::nextGaussian( int mean, int standardDeviation )
{
    return nextGaussian( standardDeviation ) + mean;
}

void Random::reset()
{
    initialize( InitialSeed );
}

unsigned int Random::initialSeed() const
{
    return InitialSeed;
}

/******************************************************************************
 * private                                                                    *
 ******************************************************************************/

unsigned int Random::nextLCPNG()
{
    // Linear Congruential Pseudo-Random Number Generator
    CurrentSeed = CurrentSeed * 0x015A4E35 + 1;
    return CurrentSeed;
}
